module drawSingle
!	!DIR$ NOOPTIMIZE
	use basicmodule
	implicit none
	
	integer		:: ir0

	contains
	
    subroutine draw_muU
		integer	:: ir,i
		real	:: m,prec,s2,v(1),u(0:q)

		prec=muprec_prior
		m=0
!		do ir=1,n
		ir=ir0
			i=clubid(ir)
			s2=om2*(1-rhos(ir)**2)*kappa2(ir)
			u=Us(:,ir)-rhos(ir)*CFs(:,i)
			m=m+sum(Siu(:,0,cindu(ir))*u)/s2
			prec=prec+Siu(0,0,cindu(ir))/s2
!		enddo
		call rnnoa(v)
		muU=m/prec+v(1)/sqrt(prec)
	end subroutine
	

	subroutine draw_rho
		integer	:: ir,i,l
		real	:: p(nrhogrid),unif(1),s2
		real	:: u(0:q)
		
!		do ir=1,n
		ir=ir0
			i=clubid(ir)
			do l=1,nrhogrid
				u=Us(:,ir)-rhogrid(l)*CFs(:,i)
				u(0)=u(0)-muU
				s2=(1-rhogrid(l)**2)*kappa2(ir)*om2
				p(l)=-.5*sum(matmul(Siu(:,:,cindu(ir)),u)*u)/s2-0.5*(q+1)*log(s2)+log(rhodist(l))
			enddo
			p=exp(p-maxval(p))
			do l=2,nrhogrid
				p(l)=p(l)+p(l-1)
			enddo
			p=p/p(nrhogrid)
			call rnun(unif)
			rhos(ir)=rhogrid(count(unif(1)>p)+1)
!		enddo
	end subroutine

	subroutine draw_CoCrho
		integer	:: i,l,k
		real	:: p(nrhogrid),unif(1),s2
		real	:: u(0:q)
		
		do i=1,1
			k=CoCid(i)
!$omp parallel do private(u,s2)			
			do l=1,nrhogrid
				u=CFs(:,i)-rhogrid(l)*CoCFs(:,k)
				s2=(1-rhogrid(l)**2)*clubs2(i)*om2
				p(l)=-.5*sum(matmul(Siu(:,:,cindCF(i)),u)*u)/s2-0.5*(q+1)*log(s2)+log(CoCrhodist(l))
			enddo
			p=exp(p-maxval(p))
			do l=2,nrhogrid
				p(l)=p(l)+p(l-1)
			enddo
			p=p/p(nrhogrid)
			call rnun(unif)
			CoCrhos(i)=rhogrid(count(unif(1)>p)+1)
		enddo
	end subroutine
	
	subroutine draw_ind(nmeas,measind,meas,measdist,meass2)
		real	:: meas(0:,:),measdist(:),meass2(:)
		integer	:: measind(:),nmeas
		
		real	:: p(nth),unif(1),u(0:q),s2
		integer	:: i,j
		
		do i=1,nmeas
			u=meas(:,i)
			s2=om2*meass2(i)
!$omp parallel do 			
			do j=1,nth
				p(j)=-.5*sum(matmul(Siu(:,:,j),u)*u)/s2+log(measdist(j)/sqrt(detSu(j)))
			enddo
			p=exp(p-maxval(p))
			do j=2,nth
				p(j)=p(j-1)+p(j)
			enddo
			p=p/p(nth)	
			call rnun(unif)
			measind(i)=count(unif(1)>p)+1
		enddo
		
	end subroutine

	subroutine draw_cindu
		real	:: meas(0:q,1)
		integer	:: ir,i

!		do ir=1,n
		ir=ir0
			i=clubid(ir)
			meas(:,1)=Us(:,ir)-rhos(ir)*CFs(:,i)
			meas(0,1)=meas(0,1)-muU
!		enddo
		call draw_ind(1,cindu(ir:ir),meas,thdist,kappa2(ir:ir)*(1-rhos(ir:ir)**2))
	end subroutine

	subroutine draw_cindCF
		real	:: meas(0:q,1)
		integer	:: i,k
		do i=1,1
			k=CoCid(i)
			meas(:,i)=CFs(:,i)-CoCrhos(i)*CoCFs(:,k)
		enddo
		call draw_ind(1,cindCF,meas,CFthdist,clubs2*(1-CoCrhos**2))
	end subroutine

	subroutine draw_cindCoC
		call draw_ind(1,cindCoC,CoCFs,CoCthdist,CoCs2)
	end subroutine

	subroutine draw_kappa2grid
		real	:: s2,u(0:q),p(ns2grid),usu,unif(1)
		integer	:: ir,i,l
		real	:: pbase(ns2grid)

		do l=1,ns2grid
			pbase(l)=-.5*(q+1)*log(s2grid(l))+log(s2dist(l))
		enddo
!		do ir=1,n
		ir=ir0
			i=clubid(ir)
			u=Us(:,ir)-rhos(ir)*CFs(:,i)
			u(0)=u(0)-muU
			s2=om2*(1-rhos(ir)**2)
			usu=sum(matmul(Siu(:,:,cindu(ir)),u)*u)
			do l=1,ns2grid
				p(l)=-.5*usu/(s2*s2grid(l))
			enddo
			p=p+pbase
			p=exp(p-maxval(p))
			do l=2,ns2grid
				p(l)=p(l)+p(l-1)
			enddo
			p=p/p(ns2grid)
			call rnun(unif)
			kappa2(ir)=s2grid(count(unif(1)>p)+1)
!		enddo
	end subroutine

	subroutine draw_clubs2grid
		real	:: s2,u(0:q),p(ns2grid),usu,unif(1)
		integer	:: i,k,l
		real	:: pbase(ns2grid)

		do l=1,ns2grid
			pbase(l)=-.5*(q+1)*log(s2grid(l))+log(CFs2dist(l))
		enddo
		do i=1,1
			k=CoCid(i)
			u=CFs(:,i)-CoCrhos(i)*CoCFs(:,k)
			s2=om2*(1-CoCrhos(i)**2)
			usu=sum(matmul(Siu(:,:,cindCF(i)),u)*u)
			do l=1,ns2grid
				p(l)=-.5*usu/(s2*s2grid(l))
			enddo
			p=p+pbase
			p=exp(p-maxval(p))
			do l=2,ns2grid
				p(l)=p(l)+p(l-1)
			enddo
			p=p/p(ns2grid)
			call rnun(unif)
			clubs2(i)=s2grid(count(unif(1)>p)+1)
		enddo
	end subroutine

	subroutine draw_CoCs2grid
		real	:: s2,u(0:q),p(ns2grid),usu,unif(1)
		integer	:: k,l
		real	:: pbase(ns2grid)

		do l=1,ns2grid
			pbase(l)=-.5*(q+1)*log(s2grid(l))+log(COCs2dist(l))
		enddo
		do k=1,1
			u=CoCFs(:,k)
			s2=om2
			usu=sum(matmul(Siu(:,:,cindCoC(k)),u)*u)
			do l=1,ns2grid
				p(l)=-.5*usu/(s2*s2grid(l))
			enddo
			p=p+pbase
			p=exp(p-maxval(p))
			do l=2,ns2grid
				p(l)=p(l)+p(l-1)
			enddo
			p=p/p(ns2grid)
			call rnun(unif)
			CoCs2(k)=s2grid(count(unif(1)>p)+1)
		enddo
	end subroutine
	
	
	subroutine draw_CF
		real	:: u(0:q),s2,Sigis(0:q,0:q,nclubs),ms(0:q,nclubs)
		integer	:: ir, i, k
		
		do i=1,1
			k=CoCid(i)
			Sigis(:,:,i)=Siu(:,:,cindCF(i))/((1-CoCrhos(i)**2)*om2*clubs2(i))
			ms(:,i)=matmul(Sigis(:,:,i),CoCrhos(i)*CoCFs(:,k))
		enddo
!		do ir=1,n
		ir=ir0
			i=clubid(ir)
			u=Us(:,ir)
			u(0)=u(0)-muU
			s2=om2*(1-rhos(ir)**2)*kappa2(ir)
			Sigis(:,:,i)=Sigis(:,:,i)+rhos(ir)**2*SiU(:,:,cindu(ir))/s2
			ms(:,i)=ms(:,i)+rhos(ir)*matmul(SiU(:,:,cindu(ir)),u)/s2
!		enddo
		do i=1,1
			call linds(Sigis(:,:,i),Sigis(:,:,i))
			call rnnoa(CFs(:,i))
			CFs(:,i)=matmul(choleski(Sigis(:,:,i)),CFs(:,i))+matmul(Sigis(:,:,i),ms(:,i))
		enddo
	end subroutine			

	subroutine draw_CoCF
		real	:: u(0:q),s2,Sigis(0:q,0:q,nCoCs),ms(0:q,nCoCs)
		integer	:: i, k
		do k=1,1
			Sigis(:,:,k)=Siu(:,:,cindCoC(k))/(om2*CoCs2(k))
			ms(:,k)=0
		enddo
		do i=1,1
			k=CoCid(i)
			u=CFs(:,i)
			s2=om2*(1-CoCrhos(i)**2)*clubs2(i)
			Sigis(:,:,k)=Sigis(:,:,k)+CoCrhos(i)**2*SiU(:,:,cindCF(i))/s2
			ms(:,k)=ms(:,k)+CoCrhos(i)*matmul(SiU(:,:,cindCF(i)),u)/s2
		enddo
		do k=1,1
			call linds(Sigis(:,:,k),Sigis(:,:,k))
			call rnnoa(CoCFs(:,k))
			CoCFs(:,k)=matmul(choleski(Sigis(:,:,k)),CoCFs(:,k))+matmul(Sigis(:,:,k),ms(:,k))
		enddo
	end subroutine			

	subroutine draw_sigma2_F
		real	:: ssum
		real	:: snu
		real	:: v(1),u(0:q),p(ns2fm),usu,unif(1)
		integer	:: l
		
		usu=sum(fm*matmul(Sifm(:,:,cindf),fm))
		do l=1,ns2fm
			p(l)=-.5*usu/s2fmtab(l)-.5*(q+1)*log(s2fmtab(l))
		enddo
		p=exp(p-maxval(p))
		do l=1,ns2fm
			p(l)=p(l)*merge(l,ns2fm+1-l,2*l<=ns2fm)
		enddo
		do l=2,ns2fm
			p(l)=p(l)+p(l-1)
		enddo
		p=p/p(ns2fm)
		call rnun(unif)
		sigma2_F(1)=s2fmtab(count(unif(1)>p)+1)
		
		u=f-fm
		u(0:1)=u(0:1)-mt_F
		ssum=sigma2_Faprior+sum(u*matmul(Sifa(:,:),u))
		snu=sigma2_Fanuprior+q+1
		call rnchi(snu,v)
		sigma2_F(2)=ssum/v(1)
	end subroutine

	
	subroutine draw_cindf
		real	:: p(nf),unif(1),u(0:q)
		integer	:: i
		
		u=fm
!		u(0:1)=u(0:1)-mt_F
!$omp parallel do 
		do i=1,nf
			p(i)=exp(-(0.5/sigma2_F(1))*sum(matmul(Sifm(:,:,i),u)*u))/sqrt(detSfm(i))
		enddo
		do i=2,nf
			p(i)=p(i-1)+p(i)
		enddo
		p=p/p(nf)
		call rnun(unif)
		cindf=count(unif(1)>p)+1
	end subroutine	

	subroutine draw_om2
		real	:: ssum,snu,v(1),s2,u(0:q)
		integer	:: ir,i,k

		ssum=om2prior
		snu=om2nuprior
!		do ir=1,n
		ir=ir0
			i=clubid(ir)
			u=Us(:,ir)-rhos(ir)*CFs(:,i)
			u(0)=u(0)-muU
			s2=(1-rhos(ir)**2)*kappa2(ir)
			ssum=ssum+sum(matmul(Siu(:,:,cindu(ir)),u)*u)/s2
			snu=snu+q+1
!		enddo
		
		do i=1,1
			k=CoCid(i)
			u=CFs(:,i)-CoCrhos(i)*CoCFs(:,k)
			s2=clubs2(i)*(1-CoCrhos(i)**2)
			ssum=ssum+sum(matmul(Siu(:,:,cindCF(i)),u)*u)/s2
			snu=snu+q+1
		enddo
		do k=1,1			
			u=CoCFs(:,k)
			s2=CoCs2(k)
			ssum=ssum+sum(matmul(Siu(:,:,cindCoC(k)),u)*u)/s2
			snu=snu+q+1
		enddo
		call rnchi(snu,v)
		om2=ssum/v(1)
	end subroutine
	
	
	subroutine draw_fcsts
		integer	:: ir,i,k
		real	:: fcstU(nh),u(0:q),s2
		
		call rnnoa(fcstf)
		fcstf=matmul(Xfcstf,mt_F)+matmul(mfcstfm(:,:,cindf),fm)+sqrt(sigma2_F(1))*matmul(cholfcstfm(:,:,cindf),fcstf)
		call rnnoa(fcstU)
		u=f-fm
		u(0:1)=u(0:1)-mt_F
		fcstf=fcstf+matmul(mfcstfa,u)+sqrt(sigma2_f(2))*matmul(cholfcstfa,fcstu)
		do k=1,1
			call rnnoa(fcstU)
			s2=om2*CoCs2(k)
			fcstCoCFs(:,k)=matmul(mfcstu(:,:,cindCoC(k)),CoCFs(:,k))+sqrt(s2)*matmul(cholfcstu(:,:,cindCoC(k)),fcstU)
		enddo
		do i=1,1
			k=CoCid(i)
			call rnnoa(fcstU)
			s2=om2*(1-CoCrhos(i)**2)*clubs2(i)
			u=CFs(:,i)-CoCrhos(i)*CoCFs(:,k)
			fcstCFs(:,i)=CoCrhos(i)*fcstCoCFs(:,k)+matmul(mfcstu(:,:,cindCF(i)),u)+sqrt(s2)*matmul(cholfcstu(:,:,cindCF(i)),fcstU)
		enddo
		
!		do ir=1,n
		ir=ir0
			i=clubid(ir)
			call rnnoa(fcstU)
			u=Us(:,ir)-rhos(ir)*CFs(:,i)
			u(0)=u(0)-muU
			s2=(1-rhos(ir)**2)*kappa2(ir)*om2
			fcstU=matmul(mfcstu(:,:,cindu(ir)),u)+sqrt(s2)*matmul(cholfcstu(:,:,cindu(ir)),fcstU)
			fcstU=fcstU+rhos(ir)*fcstCFs(:,i)
			fcstUs(:,ir)=fcstU
			fcsts(:,ir)=fcstf+fcstUs(:,ir)
!		enddo
	end subroutine
	
	subroutine draw_F
		real	:: m(0:q), u(0:q), Sigi(0:q,0:q), s2, fhat(0:q)
		integer	:: ir,i
		m=0; fhat=0
		Sigi=sigma2_f(1)*Sfm(:,:,cindf)+sigma2_f(2)*Sfa
		call linds(Sigi,Sigi)
!		do ir=1,n
		ir=ir0
			i=clubid(ir)
			s2=(1-rhos(ir)**2)*kappa2(ir)*om2
			u=Xs(:,ir)-rhos(ir)*CFs(:,i)
			u(0)=u(0)-muU
			u(0:1)=u(0:1)-mt_F
			m=m+matmul(Siu(:,:,cindu(ir)),u)/s2
			Sigi=Sigi+Siu(:,:,cindu(ir))/s2
			!if(fweights(ir)>0.0) then
			!	fhat=fhat+fweights(ir)*Xs(:,ir)
			!endif			
!		enddo
		!Sigi=Sigi+Deltainv
		!fhat(0:1)=fhat(0:1)-mt_F
		!m=m+matmul(Deltainv,fhat-Ydummy)	
		call linds(Sigi,Sigi)
		call rnnoa(f)
		f=matmul(Sigi,m)+matmul(choleski(Sigi),f)
		f(0:1)=f(0:1)+mt_F
!		do ir=1,n
		ir=ir0
			Us(:,ir)=Xs(:,ir)-F
!		enddo
	end subroutine
	
	subroutine draw_mt_F
		real	:: m(2),prec(2,2),Sigi(0:q,0:q)
		Sigi=sigma2_f(1)*Sfm(:,:,cindf)+sigma2_f(2)*Sfa
		call linds(Sigi,Sigi)
		prec=muprec_prior*eye(2)
		prec=prec+Sigi(0:1,0:1)
		m=matmul(Sigi(0:1,:),f)
		call linds(prec,prec)
		call rnnoa(mt_F)
		mt_F=matmul(prec,m)+matmul(choleski(prec),mt_f)
	end subroutine
		
	subroutine draw_fm
		real	:: u(0:q),S(0:q,0:q),Sigi(0:q,0:q),mfm(0:q,0:q)
		S=sigma2_f(1)*Sfm(:,:,cindf)
		Sigi=S+sigma2_f(2)*Sfa
		call linds(Sigi,Sigi)
		u=f
		u(0:1)=u(0:1)-mt_F
		mfm=matmul(S,Sigi)
		call rnnoa(fm)
		fm=matmul(mfm,u)+matmul(choleski(S-matmul(mfm,S)),fm)
	end subroutine

	subroutine draw_X
		integer	:: ir
		real	:: m(0:q),s2
		integer	:: i
		real	:: fhat(0:q),sVs(0:q,0:q),e(0:q),u(0:q)
		
		fhat=0;sVs=0
!		do ir=1,n
		ir=ir0
			i=clubid(ir)
			s2=(1-rhos(ir)**2)*kappa2(ir)*om2
			m=rhos(ir)*CFs(:,i)
			m(0)=m(0)+muU
			call rnnoa(u)
			u=sqrt(s2)*matmul(cholSu(:,:,cindu(ir)),u)+m
			Us(:,ir)=u-matmul(SuAA(:,:,cindu(ir),ir),u-Us(:,ir))
			!if(fweights(ir)>0.0) then
			!	fhat=fhat+fweights(ir)*Us(:,ir)
			!	sVs=sVs+fweights(ir)**2*s2*SuAAS(:,:,cindu(ir),ir)
			!endif
!		enddo
		!call rnnoa(e)
		!e=Ydummy+matmul(choleski(Delta),e)
		!fhat=matmul(invertpd(sVs+Delta),fhat-e)
!		do ir=1,n
		ir=ir0
			!if(fweights(ir)>0.0) then
			!	s2=(1-rhos(ir)**2)*kappa2(ir)*om2
			!	Us(:,ir)=Us(:,ir)-fweights(ir)*s2*matmul(SuAAS(:,:,cindu(ir),ir),fhat)
			!endif
			Xs(:,ir)=Us(:,ir)+F
!		enddo
	end subroutine

end module